AF wavefunction differential imaging dev

Exploring differential imaging, aligned vs. unaligned samples. Currently computed via effective AF wavefunctions, expanded in spherical harmonics (see Formalism below), but may not be the best way (should be able to do via matrix elements directly?).

For final phase results, see Comparison with expt section.

22/09/23

Tidy version for distro - final results only.

07/09/23

Mainly working from previous AF notebooks, see http://jake/jupyter/user/paul/doc/tree/projects-share/carlos_tross_HHG_phases_2023/notebooks/AF_wavefns_method_dev_010221_class.ipynb

  • _basic-demo version, for N2 with test alignment only. (Working with 3d-AFPAD-dev branch, requires padPlot() fix as per this commit, testing on Jake in epsdev-shared-100122 kernel.)

NOTE: currently using ePSproc base class, should update to PEMtk base.

Setup

In [1]:
# Imports
import numpy as np
import pandas as pd
import xarray as xr

# Special functions
# from scipy.special import sph_harm
import spherical_functions as sf
import quaternion

# Performance & benchmarking libraries
# from joblib import Memory
# import xyzpy as xyz
import numba as nb

# Timings with ttictoc or time
# https://github.com/hector-sab/ttictoc
# from ttictoc import TicToc
import time

# Package fns.
# For module testing, include path to module here
import sys
import os
# modPath = r'D:\code\github\ePSproc'  # Win test machine
# modPath = r'/home/femtolab/github/ePSproc/'  # Linux test machine
# sys.path.append(modPath)
import epsproc as ep
# TODO: tidy this up!
from epsproc.util import matEleSelector
from epsproc.geomFunc import geomCalc

from epsproc.classes.base import ePSbase  # Data class
from epsproc.classes.multiJob import ePSmultiJob
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
In [2]:
# Plotters
import matplotlib.pyplot as plt
from epsproc.plot import hvPlotters
hvPlotters.setPlotters()
* Set Holoviews with bokeh.
In [3]:
# Skip warnings for old XR versions
import warnings

warnings.filterwarnings('ignore', category=FutureWarning)  #message='model will result in indexing errors*')

Load data

Testing for $N_2$.

In [4]:
# Load data from modPath\data
from pathlib import Path
modPath=Path(ep.__path__[0]).parent
dataPath = os.path.join(modPath, 'data', 'photoionization')
dataFile = os.path.join(dataPath, 'n2_3sg_0.1-50.1eV_A2.inp.out')  # Set for sample N2 data for testing

data = ePSbase(fileIn = dataFile, verbose = 1)
# data = ePSbase(dataPath, verbose = 1)
data.scanFiles()

# N2O data
# dataPath = r'D:\VMs\Share\ePSshare\N2O\N20_wf'  # Test dir on Stimpy (Win machine)
# data = ePSmultiJob(dataPath, verbose = 0)
# keys = [0,1,2]  # Set for 1s datasets only (ugly!)
# data.scanFiles(keys = keys)

# # Scan data file
# dataSet = ep.readMatEle(fileIn = dataFile)
# dataXS = ep.readMatEle(fileIn = dataFile, recordType = 'CrossSection')  # XS info currently not set in NO2 sample file.
*** Job orb5 details
Key: orb5
Dir /mnt/femtobackSSHFS/DriveSyncShare/code-share/github-share/ePSproc/data/photoionization, 1 file(s).
{   'batch': 'ePS n2, batch n2_3sg_0.1-50.1eV, orbital A2',
    'event': ' N2 X-state (3sg-1)',
    'orbE': -17.34181645456815,
    'orbLabel': '3sg-1'}
In [5]:
data.jobsSummary()
*** Job orb5 details
Key: orb5
Dir /mnt/femtobackSSHFS/DriveSyncShare/code-share/github-share/ePSproc/data/photoionization, 1 file(s).
{   'batch': 'ePS n2, batch n2_3sg_0.1-50.1eV, orbital A2',
    'event': ' N2 X-state (3sg-1)',
    'orbE': -17.34181645456815,
    'orbLabel': '3sg-1'}

Matrix elements

In [6]:
# key = 'orb1' # Set for O1s
# key = 'orb2' # Set for N1s (central)
# key = 'orb3' # Set for N1s (terminal)
data.calcOpts['thres'] = 1e-4  # TODO: propagate global threshold setting! Used for calcs but not for plots at the moment.

# Plot matrix elements
# data.lmPlot(keys = key)
data.lmPlot()
Plotting data n2_3sg_0.1-50.1eV_A2.inp.out, pType=a, thres=0.01, with Seaborn
/home/paul/anaconda3/envs/epsdev-shared-100122/lib/python3.7/site-packages/xarray/core/nputils.py:215: RuntimeWarning: All-NaN slice encountered
  result = getattr(npmodule, name)(values, axis=axis, **kwargs)

Default case, unaligned

In this case the axis distribution is isotropic, and there is only a single term with $A_{Q,S}^{K}=A_{0,0}^{0}=1$.

We'll calculate all the terms in the expression for $^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})$, but skip the summations and multiplication by $Y_{l\Lambda}^{*}(\hat{k}_{M})$, thus defining effective LF/AF matrix elements.

In [7]:
# Default alignment terms, 
AKQS = ep.setADMs()
# ep.lmPlot(AKQS, xDim='t')
AKQS
Out[7]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'ADM'
  • ADM: 1
  • t: 1
  • 1
    array([[1]])
    • ADM
      (ADM)
      MultiIndex
      (K, Q, S)
      array([(0, 0, 0)], dtype=object)
    • K
      (ADM)
      int64
      0
      array([0])
    • Q
      (ADM)
      int64
      0
      array([0])
    • S
      (ADM)
      int64
      0
      array([0])
    • t
      (t)
      int64
      0
      units :
      Index
      array([0])
  • dataType :
    ADM
    long_name :
    Axis distribution moments
    units :
    arb
In [8]:
# Run default af numeric code
# This will default to Type='L', it=1 and an isotropic axis distribution.
data.afpadNumeric()
/home/paul/anaconda3/envs/epsdev-shared-100122/lib/python3.7/site-packages/xarray/core/nputils.py:215: RuntimeWarning: All-NaN slice encountered
  result = getattr(npmodule, name)(values, axis=axis, **kwargs)
Set AF expansion parameters TXaf and DeltaKQS for orb5
In [9]:
# Set as ref data
key = 'orb5'
for item in ['TXaf', 'DeltaKQS']:
    data.data[key][item + 'Unaligned']= data.data[key][item]

With alignment

Here we can look at the dependence of the channel coupling parameters $^{AF}\Delta_{l,m}(K,Q,S)$, as well as the final results for a given axis distribution...

Channel couplings for various $K$

In this case, assume a $z$-polarized ionizing field, and a cylindrically symmetric distribution. To look at the channel couplings, set all terms $A_{Q,S}^{K}=A_{0,0}^{K}=1$ for even-$K$ up to $K_{max}=8$.

In [10]:
# # Set ADMs - TEST ONLY
# KQSlist = [[K,0,0,1] for K in range(0,9,2)]
# AKQS = ep.setADMs(ADMs = KQSlist)  # TODO: routine to automate this type of setting!

# # daPlot, daPlotpd, legendList, gFig = ep.lmPlot(AKQS, xDim = 't', pType = 'r', squeeze = False)  # Note squeeze = False required for 1D case (should add this to code!)
# # daPlotpd
# AKQS #.coords

# SET ADMs from file, see https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_basic_demo_030621-full_010922.html
# Load time-dependent ADMs for N2 case
# Adapted from ePSproc_AFBLM_testing_010519_300719.m

from scipy.io import loadmat
ADMdataFile = os.path.join(modPath, 'data', 'alignment', 'N2_ADM_VM_290816.mat')
ADMs = loadmat(ADMdataFile)

# Set tOffset for calcs, 3.76ps!!!
# This is because this is 2-pulse case, and will set t=0 to 2nd pulse (and matches defn. in N2 experimental paper)
tOffset = -3.76
ADMs['time'] = ADMs['time'] + tOffset

AKQS2pulse = ep.setADMs(ADMs = ADMs['ADM'], t=ADMs['time'].squeeze(), KQSLabels = ADMs['ADMlist'], addS = True)

# Subselect
# AKQS2pulseSub = AKQS2pulse.sel(t=slice(4, 5, 4))
AKQS2pulseSub = AKQS2pulse.sel(t=slice(4, 6.5, 4))
# AKQS2pulseSub

# PEMtk class has this directly:
# data.setADMs(ADMs = ADMs['ADM'], t=ADMs['time'].squeeze(), KQSLabels = ADMs['ADMlist'], addS = True)
# data.data['ADM']['ADM']
# data.selOpts['ADM'] = {}   #{'thres': 0.01, 'inds': {'Type':'L', 'Eke':1.1}}
# data.setSubset(dataKey = 'ADM', dataType = 'ADM', sliceParams = {'t':[4, 5, 4]})
# data.ADMplot(keys = 'subset')

Aligned frame angular distributions

Finally, the observable angular patterns can be investigated, by looking at $|^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})|^2$.

NOTE data types for plots:

  • 'r': real
  • 'i': imaginary
  • 'a': abs values
  • 'phaseUW': phase, unwrapped (by energy)
In [11]:
# Run AF numeric code with unity AKQS parameters
data.afpadNumeric(AKQS = AKQS2pulseSub)
Set AF expansion parameters TXaf and DeltaKQS for orb5
/home/paul/anaconda3/envs/epsdev-shared-100122/lib/python3.7/site-packages/xarray/core/nputils.py:215: RuntimeWarning: All-NaN slice encountered
  result = getattr(npmodule, name)(values, axis=axis, **kwargs)
In [12]:
# Plot panels with padPlot
Erange = [1,40,1]

# This currently doesn't support the native dims, so reset these first and stack to standard LM representation (includes sum over MF terms)
for key in data.data.keys():
    data.data[key]['TXafSum'] = data.data[key]['TXaf'].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
    data.data[key]['TXafSum'].attrs = data.data[key]['TXaf'].attrs
    

data.padPlot(Etype = 'Eke', Erange = Erange, dataType='TXafSum', pType='a2', 
             backend='mpl', sumDims=['Sym'], facetDims=['t', 'Eke'], selDims={'mu0':0}, 
             reducePhi='sum', pStyle='grid', returnFlag=True)  #, squeeze=True) 

# TO TEST: should be able to add  col_wrap=3 for Xarray plotter?
# See https://docs.xarray.dev/en/latest/generated/xarray.DataArray.plot.html
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafSum], facetDims=['Eke', 't'], pType=a2 with backend=mpl.
Grid plot: 3sg-1, dataType: TXafSum, plotType: a2
Set plot to self.data['orb5']['plots']['TXafSum']['grid']

Differential case

In [13]:
# Set data
for key in data.data.keys():
    for dataType in ['TXaf','TXafUnaligned']: 
        sumData = dataType + 'Sum'
        data.data[key][sumData] = data.data[key][dataType].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
        data.data[key][sumData].attrs = data.data[key][dataType].attrs
        
        # Set distributions
        for pType in ['r','i','a','phaseUW']:
            
            if 't' in data.data[key][sumData].squeeze().dims:
                facetDims=['t', 'Eke']
            else:
                facetDims=['Eke']
                
            
            data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType=pType, 
                     backend='mpl', sumDims=['Sym'], facetDims=facetDims, selDims={'mu0':0}, 
                     reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False);  #, squeeze=True) 
            

            # Set plots by type (currently NOT done above!)
#             if not data.data[key]['plots'][sumData]:
#                 data.data[key]['plots'][sumData][pType] = {}
                
#             print(f"Setting data.data[{key}]['plots'][{sumData}][{pType}]")
#             data.data[key]['plots'][sumData][pType]['pData'] = data.data[key]['plots'][sumData]['pData'].copy()
#             data.data[key]['plots'][sumData][pType]['grid'] = data.data[key]['plots'][sumData]['grid'].copy()
    
            # Dataout
            outData = sumData + 'Out'
            if not outData in data.data[key].keys():
                data.data[key][outData] = {}
            
            print(f"Setting data.data[{key}][{outData}][{pType}]")
            data.data[key][outData][pType] = {'pData': data.data[key]['plots'][sumData]['pData'].copy()}
            
            # Subtract
        
    
#     data.data[key]['TXafUnalignedSum'] = data.data[key]['TXaf'].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafSum], facetDims=['Eke', 't'], pType=r with backend=mpl.
Grid plot: 3sg-1, dataType: TXafSum, plotType: r
Set plot to self.data['orb5']['plots']['TXafSum']['grid']
Setting data.data[orb5][TXafSumOut][r]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafSum], facetDims=['Eke', 't'], pType=i with backend=mpl.
Grid plot: 3sg-1, dataType: TXafSum, plotType: i
Set plot to self.data['orb5']['plots']['TXafSum']['grid']
Setting data.data[orb5][TXafSumOut][i]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafSum], facetDims=['Eke', 't'], pType=a with backend=mpl.
Grid plot: 3sg-1, dataType: TXafSum, plotType: a
Set plot to self.data['orb5']['plots']['TXafSum']['grid']
Setting data.data[orb5][TXafSumOut][a]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafSum], facetDims=['Eke', 't'], pType=phaseUW with backend=mpl.
Grid plot: 3sg-1, dataType: TXafSum, plotType: phaseUW
Set plot to self.data['orb5']['plots']['TXafSum']['grid']
Setting data.data[orb5][TXafSumOut][phaseUW]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=r with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: r
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
Setting data.data[orb5][TXafUnalignedSumOut][r]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=i with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: i
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
Setting data.data[orb5][TXafUnalignedSumOut][i]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=a with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: a
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
Setting data.data[orb5][TXafUnalignedSumOut][a]
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=phaseUW with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: phaseUW
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
Setting data.data[orb5][TXafUnalignedSumOut][phaseUW]
In [14]:
# # Unaligned ref:
data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType='a', 
                     backend='mpl', sumDims=['Sym'], facetDims=['Eke'], selDims={'mu0':0}, 
                     reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False);
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=a with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: a
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
In [15]:
data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType='phaseUW', 
                     backend='mpl', sumDims=['Sym'], facetDims=['Eke'], selDims={'mu0':0}, 
                     reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False);
*** sphPlot dataType = afpad not recognised, trying anyway.
Using default sph betas.
Summing over dims: {'Sym'}
Found additional dims {'t', 'E'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb5][TXafUnalignedSum], facetDims=['Eke', None], pType=phaseUW with backend=mpl.
Grid plot: 3sg-1, dataType: TXafUnalignedSum, plotType: phaseUW
Set plot to self.data['orb5']['plots']['TXafUnalignedSum']['grid']
In [38]:
# Subtract distributions
# Note this will also define the sign of the phases
subData = 'sub'
data.data[key][subData] = {}
for pType in ['r','i','a','phaseUW']:
    data.data[key][subData][pType] = {}
    data.data[key]['sub'][pType]['pData'] = data.data[key]['TXafUnalignedSumOut'][pType]['pData'] - \
                                            data.data[key]['TXafSumOut'][pType]['pData']
                                            

(E,Theta) differential plots

These should be respresentative of measurements...

NOTE data types:

  • 'r': real
  • 'i': imaginary
  • 'a': abs values
  • 'phaseUW': phase, unwrapped (by energy)
In [39]:
# PLOTS

pType = 'r'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
Out[39]:
<xarray.plot.facetgrid.FacetGrid at 0x7fcaa1632c10>
In [40]:
# PLOTS

pType = 'i'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
Out[40]:
<xarray.plot.facetgrid.FacetGrid at 0x7fcaa08e60d0>
In [41]:
# PLOTS

pType = 'a'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
Out[41]:
<xarray.plot.facetgrid.FacetGrid at 0x7fcaa15d2490>
In [42]:
# PLOTS

pType = 'phaseUW'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
Out[42]:
<xarray.plot.facetgrid.FacetGrid at 0x7fc7ef105d50>

E-only diff plots - Eke scale

Summed over angular coord...

TODO: should pass full complex form here before sum!

Current case maybe incorrect - 10^-13 for r/i components? (Two anti-phase theta components cancel!)

In [43]:
# Sum over theta...
pType = 'r'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke')  #,col='t', col_wrap=3)

plt.title(f"Data type: {pType}")
Out[43]:
Text(0.5, 1.0, 'Data type: r')
In [44]:
# Sum over theta...
pType = 'i'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke')  #,col='t', col_wrap=3)

plt.title(f"Data type: {pType}")
Out[44]:
Text(0.5, 1.0, 'Data type: i')
In [45]:
# Sum over theta...
pType = 'a'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke')  #,col='t', col_wrap=3)

plt.title(f"Data type: {pType}")
Out[45]:
Text(0.5, 1.0, 'Data type: a')
In [46]:
# Sum over theta...
pType = 'phaseUW'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke')  #,col='t', col_wrap=3)

plt.title(f"Data type: {pType}")
Out[46]:
Text(0.5, 1.0, 'Data type: phaseUW')

THINK this should roughly compare to the results in Fig. 8.7 of Jan Tross' thesis - need to convert to harmonics scale, but he sees large feature around 4 - 5ps, flips sign between H9 and H13.

E-only diff plots - Harmonics scale

Summed over angular coord...

TODO: should pass full complex form here before sum!

Current case maybe incorrect - 10^-13 for r/i components? (Two anti-phase theta components cancel!)

In [47]:
data.data['orb5']['XS'].Ehv
Out[47]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'Ehv'
  • Eke: 51
  • 17.4 18.4 19.4 20.4 21.4 22.4 23.4 ... 62.4 63.4 64.4 65.4 66.4 67.4
    array([17.4, 18.4, 19.4, 20.4, 21.4, 22.4, 23.4, 24.4, 25.4, 26.4, 27.4,
           28.4, 29.4, 30.4, 31.4, 32.4, 33.4, 34.4, 35.4, 36.4, 37.4, 38.4,
           39.4, 40.4, 41.4, 42.4, 43.4, 44.4, 45.4, 46.4, 47.4, 48.4, 49.4,
           50.4, 51.4, 52.4, 53.4, 54.4, 55.4, 56.4, 57.4, 58.4, 59.4, 60.4,
           61.4, 62.4, 63.4, 64.4, 65.4, 66.4, 67.4])
    • Ehv
      (Eke)
      float64
      17.4 18.4 19.4 ... 65.4 66.4 67.4
      array([17.4, 18.4, 19.4, 20.4, 21.4, 22.4, 23.4, 24.4, 25.4, 26.4, 27.4,
             28.4, 29.4, 30.4, 31.4, 32.4, 33.4, 34.4, 35.4, 36.4, 37.4, 38.4,
             39.4, 40.4, 41.4, 42.4, 43.4, 44.4, 45.4, 46.4, 47.4, 48.4, 49.4,
             50.4, 51.4, 52.4, 53.4, 54.4, 55.4, 56.4, 57.4, 58.4, 59.4, 60.4,
             61.4, 62.4, 63.4, 64.4, 65.4, 66.4, 67.4])
    • Eke
      (Eke)
      float64
      0.1 1.1 2.1 3.1 ... 48.1 49.1 50.1
      array([ 0.1,  1.1,  2.1,  3.1,  4.1,  5.1,  6.1,  7.1,  8.1,  9.1, 10.1, 11.1,
             12.1, 13.1, 14.1, 15.1, 16.1, 17.1, 18.1, 19.1, 20.1, 21.1, 22.1, 23.1,
             24.1, 25.1, 26.1, 27.1, 28.1, 29.1, 30.1, 31.1, 32.1, 33.1, 34.1, 35.1,
             36.1, 37.1, 38.1, 39.1, 40.1, 41.1, 42.1, 43.1, 44.1, 45.1, 46.1, 47.1,
             48.1, 49.1, 50.1])
In [48]:
# TODO: util function for this!
# See ep.util.conversion.conv_ev_nm

#     # Define constants from scipy.constants
#     h = scipy.constants.h
#     c = scipy.constants.c
#     evJ = scipy.constants.physical_constants['electron volt-joule relationship'][0]

#     # Define output units - wavelength in m
#     waveConv = 1e-9
#     dataOut = (h * c)/(data * evJ)/waveConv

import scipy.constants
def setHHG(dataIn, w=785, Etype='hv', Eshift=0, dp=2):
    """
    Set harmonics based on driving field photon energy.
    
    Assumes Xarray input, and sets according to dataIn.Ehv/photonE.
    
    Set a shift if required (since dataIn.Ehv may be off from real values).
    
    dp is used to round values in output.
    
    """
    
#     h = scipy.constants.h
#     c = scipy.constants.c
#     evJ = scipy.constants.physical_constants['electron volt-joule relationship'][0]
#     waveConv = 1e-9
    
#     # Photon energy
#     photonE = (h * c)/(w * evJ)/waveConv

    photonE = ep.util.conversion.conv_ev_nm(w)
    
    print(f"Converting to harmonic E scale with \Omega={w} nm, photon E={np.round(photonE,3)} eV from Ehv, Eshift={Eshift}")
    
    # Convert data scale to harmonics of specified E
    dataIn['HH'] = ((dataIn.Ehv+Eshift)/photonE).round(dp)
    
    dataIn.attrs['HHGconv'] = {'HH':dataIn['HH'].data.copy(),
                               'photonE':photonE,
                               '\omega':w}
    
    print(f"Photon Erange [{dataIn.Ehv.data[0]}:{dataIn.Ehv.data[-1]}] set to harmonics \
            [{np.round(dataIn.HH.data[0],2)}:{np.round(dataIn.HH.data[-1],2)}]")
    
#     return testCoords
    
# testCoords = data.data['orb5']['XS'].Ehv.copy()

# setHHG(testCoords)
# data.data['orb5']['XS'].Ehv

# Set HHG axis in data.
# Includes Eshift to match 1st IP (14.5 eV), vs. 17.5 in orb5 results (but may be wrong orb!!!).
# TODO - check data attrs, FegeEng vs. orbIP
# setHHG(data.data['orb5']['XS'], Eshift=-3)

for pType in data.data[key][subData].keys():
    setHHG(data.data[key][subData][pType]['pData'], Eshift=-4.2)
Converting to harmonic E scale with \Omega=785 nm, photon E=1.579 eV from Ehv, Eshift=-4.2
Photon Erange [18.4:56.4] set to harmonics             [8.99:33.05]
Converting to harmonic E scale with \Omega=785 nm, photon E=1.579 eV from Ehv, Eshift=-4.2
Photon Erange [18.4:56.4] set to harmonics             [8.99:33.05]
Converting to harmonic E scale with \Omega=785 nm, photon E=1.579 eV from Ehv, Eshift=-4.2
Photon Erange [18.4:56.4] set to harmonics             [8.99:33.05]
Converting to harmonic E scale with \Omega=785 nm, photon E=1.579 eV from Ehv, Eshift=-4.2
Photon Erange [18.4:56.4] set to harmonics             [8.99:33.05]
In [49]:
# Sum over theta... REPLOT ON HH SCALE
pType = 'phaseUW'
# pType = 'a'

# data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sum('Theta').real.plot(y='HH')
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='HH')

plt.title(f"Data type: {pType}")
Out[49]:
Text(0.5, 1.0, 'Data type: phaseUW')
In [50]:
# Pull discrete harmonics...
HHaxis = np.arange(9,21,2)
data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sel(HH=HHaxis, method="nearest") \
    .sum('Theta').real.plot(x='t')

plt.title(f"Data type: {pType}")
Out[50]:
Text(0.5, 1.0, 'Data type: phaseUW')

Comparison with expt

In [51]:
# Line plot
data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sel(HH=HHaxis, method="nearest") \
    .sum('Theta').real.plot.line(x='t')

plt.title(f"Data type: {pType}")
Out[51]:
Text(0.5, 1.0, 'Data type: phaseUW')

Compare phases to Jan's thesis, Fig. 8.7.

  • Similar behaviour.
  • Switch in sign of phases between H9 and H13.
  • t-axis may be slightly different in these cases (calcs herein are for 2-pulse alignment case).

Looks reasonable, given that the calcs herein use slightly different alignment parameters, and are for the HOMO only.

In [58]:
from IPython.display import Image
from IPython.core.display import HTML 
# Image(url= "N2_phases_tross_thesis_fig8.7_crop.png", width=500, height=100)
Image(url= "N2_phases_tross_thesis_fig8.7.png", width=800, height=100)
Out[58]:

Formalism

(TO CHECK - conj D terms...?)

The usual ePS eqns. for the MF results express the MF wavefunction as an expansion in the radial matrix elements plus some angular terms:

\begin{equation} T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\sum_{l,m,\mu}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})D_{\mu,\mu_{0}}^{1}(R_{\hat{n}}) \end{equation}

The MFPAD is then determined from the square of this function. Note here that the polarisation term is usually expressed in the LF, then rotated into the MF.

The MFPAD can also be expressed in terms of $\beta_{LM}$ parameters, derived directly from the analytical square of the wavefunction. In this case the various angular momenta are coupled directly, and the observable is expressed as a set of $\beta_{LM}$s. This is usually preferable for computation of observables, however, since these are derived from the coherent square of the wavefunction, there is no phase information, so for coupling into further calculations one may prefer to work with the effective matrix elements or wavefunctions directly.

For the LF the same considerations apply, although one usually works with the $\beta_{LM}$s. For the LF the wavefunctions can be written as:

\begin{equation} ^{LF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L},R_{\hat{n}})=\sum_{l,m,\mu,\Lambda}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)D_{m\Lambda}^{l*}(\Omega)Y_{l\Lambda}^{*}(\hat{k}_{L})D_{\mu,\mu_{0}}^{1*}(\Omega) \end{equation}

where the notation has been simplified (following [1]), and an additional term $D_{m\Lambda}^{l}(\Omega)$ (where $\Omega$ denotes Euler angles $(\phi,\theta,\chi)$) rotates the MF projection terms $m$ into the LF, and all Euler angles are integrated over (i.e. all molecular alignments/MF polarization orientations). Note this is usually written with the opposite labels, i.e. $m$ is the LF projection term and $\lambda$ the MF term, but we'll keep the ePS notation for consistency here.

For the AF, the distribution is essentially a convolution (assuming that the rotational wavefunction is separable). (See the AFBLM pages for more.)

ROUGHLY:

\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L},R_{\hat{n}})=\intop d\Omega\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)D_{m\Lambda}^{l*}(\Omega)Y_{l\Lambda}^{*}(\hat{k}_{M})D_{\mu,\mu_{0}}^{1*}(\Omega)A_{Q,S}^{K}D_{Q,S}^{K*}(\Omega) \end{equation}

Where $A_{Q,S}^{K}D_{Q,S}^{K}(\Omega)$ describes the axis distribution function in the LF Euler angles.

Here there is a triple product in $D$ - convert to non-conj terms (unnecessary?) then apply eqn. 3.118 in Zare:

\begin{eqnarray} \intop d\Omega D_{m\Lambda}^{l*}(\Omega)D_{\mu,\mu_{0}}^{1*}(\Omega)D_{Q,S}^{K*}(\Omega) & = & \intop d\Omega(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}D_{-m,-\Lambda}^{l}(\Omega)D_{-\mu,-\mu_{0}}^{1}(\Omega)D_{-Q,-S}^{K}(\Omega)\\ & = & 8\pi^{2}(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right) \end{eqnarray}

Hence:

\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})=8\pi^{2}\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}A_{Q,S}^{K}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right)Y_{l\Lambda}^{*}(\hat{k}_{M}) \end{equation}

For the geometric computations, define this similar to existing AF $\Delta_{L,M}(K,Q,S)$ term (for AF $\beta_{LM}$ calculations):

\begin{equation} ^{AF}\Delta_{l,m}(K,Q,S)=(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right) \end{equation}

Hence:

\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})=8\pi^{2}\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}A_{Q,S}^{K}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E){}^{AF}\Delta_{l,m}(K,Q,S)Y_{l\Lambda}^{*}(\hat{k}_{M}) \end{equation}

Where $^{AF}\Delta_{l,m}(K,Q,S)$ is a channel coupling parameter.

[1] Underwood, Jonathan G., and Katharine L. Reid. “Time-Resolved Photoelectron Angular Distributions as a Probe of Intramolecular Dynamics: Connecting the Molecular Frame and the Laboratory Frame.” The Journal of Chemical Physics 113, no. 3 (2000): 1067. https://doi.org/10.1063/1.481918.

NEXT:

  • Check result above more carefully, may have messed up some phases (likely). Should be non-conj term in $\mu,\mu_0$?
  • Test numerically - should be able to compare numerical vs. beta results, similarly to LF test cases

Versions

In [31]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
Out[31]:
Fri Sep 22 10:59:22 2023 EDT
OS Linux CPU(s) 64 Machine x86_64
Architecture 64bit Environment Jupyter
Python 3.7.10 | packaged by conda-forge | (default, Feb 19 2021, 16:07:37) [GCC 9.3.0]
epsproc 1.3.2-dev xarray 0.15.0 jupyter Version unknown
numpy 1.21.6 scipy 1.6.1 IPython 7.21.0
matplotlib 3.3.4 scooby 0.5.6
In [32]:
# Check current Git commit for local ePSproc version
from pathlib import Path
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
* 3d-AFPAD-dev
  dev
  master
c1eedf10631cb17fe33421427c195a2391369c37
In [33]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
c1eedf10631cb17fe33421427c195a2391369c37	refs/heads/3d-AFPAD-dev
897d73392a7b32ffba4ca6b6b4755c61e7c1c8d7	refs/heads/dependabot/pip/notes/envs/envs-versioned/certifi-2022.12.7
457f8cd85d89bd6474296b6c01e5165a4a7ce7fc	refs/heads/dependabot/pip/notes/envs/envs-versioned/cryptography-39.0.1
2855573d0f088b45d19acf2fd9a71eeb7af0a29b	refs/heads/dependabot/pip/notes/envs/envs-versioned/ipython-8.10.0
92c661789a7d2927f2b53d7266f57de70b3834fa	refs/heads/dependabot/pip/notes/envs/envs-versioned/mistune-2.0.3
fe1e9540c7b91fe571f60562acd31d8e489d491e	refs/heads/dependabot/pip/notes/envs/envs-versioned/nbconvert-6.5.1
70b80a1e3a54de91c2bfe3b6be82d611fcfd5f43	refs/heads/dependabot/pip/notes/envs/envs-versioned/pillow-9.3.0
92fc79b09aafedadcb645f88bb7ed771c96d5b52	refs/heads/dependabot/pip/notes/envs/envs-versioned/setuptools-65.5.1
fa33ed8d63a5c4a4043cc4c261059cc09e4c2bf7	refs/heads/dependabot/pip/notes/envs/envs-versioned/wheel-0.38.1
41cdfe43750e08c510f98b05e024a9c62da42771	refs/heads/dependabot/pip/setuptools-65.5.1
7e4270370d66df44c334675ac487c87d702408da	refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master
baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7	refs/heads/testDev
In [ ]: